/*
 * graph.c
 *
 * this module provides functions to manipulate the adjacency matrix of a graph.
 *
 *  Created on: Aug 16, 2011
 *      Author: correa
 */

#include <cstdio>
#include <cstdlib>
#include "heap.h"
#include "Graph.h"
#include <gmp.h>


#define declareListForAssignment(r,k) 	int _i_; \
		int _next_[r+1]; \
		for (_i_ = 0; _i_ <= r; _i_++) _next_[_i_] = r; \
		int _covhead_[k+1]; \
		for (_i_ = 0; _i_ <= k; _i_++) _covhead_[_i_] = r

#define callAssignClique(g, k, cover, i, coversize)	assignCliqueByMatrix(g, k, cover, i, &coversize, _next_, _covhead_)

int *satdeg;

template <class W>
static inline int assignByList(const W * d, int ** ng, const int k, int * cover, const int i, int * coversize, int * clique) {
	int j;
	// current clique starts with size 1
	cover[i] = 0;

	// find clique to include vertex as the smallest
	// not appearing in the neighborhood of i
	memset(clique, 0, ((*coversize)+1)*sizeof(int));
	for (j = 0; j < d[i]; j++)
		clique[cover[ng[i][j]]] = 1;

	int c = 1;
	while (c <= *coversize && clique[c])
		c++;
	if (c <= k) {
		if (c > *coversize)
			(*coversize)++;
		cover[i] = c;
	}

	return c;
}

template <class T>
static inline int assignCliqueByMatrix(T * const g, const int k, int * cover, const int i, int * coversize, int * next, int * covhead) {
	int c = 0, r = g->nverts();
	int cc = covhead[c];
	int antcc = r;
	while (!cover[i] && c < *coversize) {
		// search for cover to add
		for (; cc != covhead[c+1] && g->hasEdge(i, cc); antcc = cc, cc = next[cc]);
		//			printf("has edge (%d,%d)\n",g->vertex(i), g->vertex(cc));
		if (cc == covhead[c+1])
			cover[i] = c + 1;
		else {
			c++;
			for (; cc != covhead[c]; antcc = cc, cc = next[cc]);
		}
	}
	if (c == *coversize && *coversize < k) {
		(*coversize)++;
		covhead[c] = i;
		cover[i] = c + 1;
	}
	// include i in the tail of cover list
	next[antcc] = i;
	next[i] = cc;

	return c + 1;
}


static inline int max_satdeg(const void * a, const void * b) {
	if(satdeg[*((int*)a)] > satdeg[*((int*)b)]) return 1;
	else {
		if(satdeg[*((int*)a)] < satdeg[*((int*)b)]) return -1;
		else return 0;
	}
}


// sort R according to a clique cover or coloring of the subgraph induced by R,
// depending whether d and ng represents neighborhoods or antineighborhoods.
// index in R is used as last break tie condition
// Assuming they represent antineighborhoods, clique cover is determined using dsatur
// 4 vectors of r integers are defined in the C function stack
// it returns the size of the clique cover and sets the cover array such that:
// for i in 0, ..., r-1, the entry cover[i] contains the size of a clique
// cover of the subgraph induced by all vertices of rank less than or equal to cover[i]
// the number of cliques in the cover is returned. Cliques returned in the cover array are numbered from 1 to number of cliques in the
template <class W>
static inline int dsaturSort_impl(Graph * const g, Weight<W> * w, const int * d, int ** ng, const int k, int * rank) {
	int r = g->nverts();
	int next[r+1];
	int prev[r+1];
	W satdeg[r];
	int satdegheap[r];
	unsigned int nmemb = 1;
	W colorContr[r][r];

	// array, indexed by the cliques, used to mark all cliques in the neighborhood of imaxantideg
	int clique[k+1];
	int nclique = 0;

	int i, j, l,aux;

	// initially, all cliques are empty, all saturation degrees are 0 and
	// its doubly linked list contains all vertices
	for (i = 0; i < r; i++) {
		satdegheap[i]=r;
		for (j = 0; j < r; j++) {
			colorContr[i][j]=(W) 0;
		}
	}

	memset(rank, 0, r*sizeof(int));
	memset(satdeg, 0, r*sizeof(W));
	for (i = 0; i < r; i++) {
		next[i] = i+1;
		prev[i] = i-1;
		satdegheap[i] = r;
	}
	satdegheap[0]=0;
	prev[0] = r;

	int ncovered;
	int imaxantideg;
	int found = 0;
	for (ncovered = 0; ncovered < r; ncovered++) {
		// find vertex with max antideg among those with max satdeg
		i = satdegheap[0];
		imaxantideg = i;
		while (next[i] < r) {
			i = next[i];
			if (d[i] > d[imaxantideg])
				imaxantideg = i;
			else {
				 if ( (d[i] == d[imaxantideg]) && (i<imaxantideg) ) imaxantideg = i;
			}
		}

		// remove vertex form its doubly linked list
		if (prev[imaxantideg] == r) {
			satdegheap[0] = next[imaxantideg];
			if (ncovered < r - 1)
				if(satdegheap[0] == r) {
					aux=satdegheap[nmemb-1];
					satdegheap[nmemb-1]=satdegheap[0];
					satdegheap[0]=aux;
					nmemb--;
					heapdown(0, &satdegheap[nmemb], satdegheap, nmemb, sizeof(int), max_satdeg);
				}
		}
		else
			next[prev[imaxantideg]] = next[imaxantideg];
		prev[next[imaxantideg]] = prev[imaxantideg];

		// find clique to include vertex as the smallest
		// not appearing in the neighborhood of imaxantideg
		int c = assignByList(d, ng, k, rank, imaxantideg, &nclique, clique);

		// update satdeglist
		if (c <= k)
			for (i = 0; i < d[imaxantideg]; i++) {
				int ii = ng[imaxantideg][i];
				if (!rank[ii]) {
					for (j = 0; j < d[ii] && (ng[ii][j] == imaxantideg || rank[ng[ii][j]] != c); j++);
					if (j == d[ii]) {
						if(w->weight(ii,imaxantideg) > colorContr[ii][c]) {
							// remove i from list
							if (prev[i] == r) {
								for(l=0;l<nmemb;l++) {
									if( satdeg[satdegheap[l]] == satdeg[i] ) break;
								}
								if(next[i] == r) {
									aux=satdegheap[nmemb-1];
									satdegheap[nmemb-1]=satdegheap[l];
									satdegheap[l]=aux;
									nmemb--;
									if (satdeg[satdegheap[l]] > satdeg[satdegheap[l/2]]) heapup(satdegheap, l, sizeof(int), max_satdeg);
									else heapdown(l, &satdegheap[nmemb], satdegheap, nmemb, sizeof(int), max_satdeg);
								}
								else satdegheap[l] = next[i];
							} else next[prev[i]] = next[i];
							prev[next[i]] = prev[i];

							// calculate the new sat degree

							satdeg[i]+= (w->weight(i,j) - colorContr[i][c]);
							colorContr[i][c] = w->weight(i,j);
							// insert i in the list with new sat deg

							found = 0;
							for(l=0;l<nmemb;l++) {
								if(satdeg[satdegheap[l]] == satdeg[i]) {
									next[i] = satdegheap[l];
									prev[i] = r;
									prev[satdegheap[l]] = i;
									satdegheap[l] = i;
									found = 1;
									break;
								}
							}

							if(!found) {
								heapoffer(&i, satdegheap, nmemb, sizeof(int), max_satdeg);
								next[i] = r;
								prev[i] = r;
							}
						}
					}
				}
			}
	}

	return nclique;
}

template <>
int WGraphs<int>::dsaturSort(Graph * const g, Weight<int> * w, const int * d, int ** ng, const int k, int * rank) {
	return dsaturSort_impl<int>(g,w, d, ng, k, rank);
}

template <>
int WGraphs<double>::dsaturSort(Graph * const g, Weight<double> * w, const int * d, int ** ng, const int k, int * rank) {
	return dsaturSort_impl<double>(g,w, d, ng, k, rank);
}

// sort R according to a clique cover of the subgraph induced by R,
// using the symmetric adjacency matrix
// clique cover is determined using dsatur
// 5 vectors of r integers are defined in the C function stack
// it returns the size of the clique cover and sets the cover array such that:
// for i in 0, ..., r-1, the entry cover[i] contains the size of a clique
// cover of the subgraph induced by all vertices of rank less than or equal to cover[i]
// the number of cliques in the cover is returned. Cliques returned in the cover array are numbered from 1 to number of cliques in the
template <class W>
static inline int dsaturCoverSort_impl(Graph * const g, Weight<W> * w, const int k, int * rank) {
	const int r = g->nverts();
	declareListForAssignment(r, k);

	int antideg[r];
	int next[r+1];
	int prev[r+1];
	W satdeg[r];
	int nclique = 0;
	int satdegheap[r];
	unsigned int nmemb;

	W ** colorContr = new W*[r];

	int i, j, l, aux;

	// initially, all cliques are empty, all saturation degrees are 0 and
	// its doubly linked list contains all vertices
	memset(rank, 0, r*sizeof(int));
	memset(antideg, 0, r*sizeof(int));
	for (i = 0; i < r; i++) {
		satdeg[i] = (W) 0;
		colorContr[i] = new W[r];
		for (j = 0; j < r; j++) {
			colorContr[i][j]=(W) 0;
		}
	}

	for (i = 0; i < r; i++) {
		next[i] = i+1;
		prev[i] = i-1;
		satdegheap[i] = r;
		for (j = i + 1; j < r; j++)
			if (!g->hasEdge(i, j)) {
				antideg[i]++;
				antideg[j]++;
			}
	}
	satdegheap[0]=0;
	nmemb = 1;
	prev[0] = r;

	int ncovered;
	int imaxantideg;
	int found = 0;
	for (ncovered = 0; ncovered < r; ncovered++) {
		// find vertex with max antideg among those with max satdeg
		i = satdegheap[0];
		imaxantideg = i;
		while (next[i] < r) {
			i = next[i];
			if (antideg[i] > antideg[imaxantideg])
				imaxantideg = i;
			else {
				  if ( (antideg[i] == antideg[imaxantideg]) && (i<imaxantideg) ) imaxantideg = i;
			}
		}

		// remove vertex from its doubly linked list
		if (prev[imaxantideg] == r) {
			satdegheap[0] = next[imaxantideg];
			if (ncovered < r - 1 && satdegheap[0] == r)
				if (--nmemb > 0) {
					satdegheap[0] = satdegheap[nmemb];
					heapdown(0, &satdegheap[nmemb], satdegheap, nmemb, sizeof(int), max_satdeg);
				}
		}
		else
			next[prev[imaxantideg]] = next[imaxantideg];
		prev[next[imaxantideg]] = prev[imaxantideg];

		int c = callAssignClique(g, k, rank, imaxantideg, nclique);
		//		printf("clique of %d is %d\n",g->vertex(imaxantideg),c);

		// update satdeglist
		if (c <= k)
			for (i = 0; i < r; i++)
				if (!g->hasEdge(imaxantideg, i) && !rank[i]) {
					for (j = 0; j < r && (g->hasEdge(i, j) || j == imaxantideg || rank[j] != c); j++);
					if (j == r) {
						if(w->weight(i,imaxantideg) > colorContr[i][c]) {

							// remove i from list
							if (prev[i] == r) {
								for(l=0;l<nmemb;l++) {
									if( satdeg[satdegheap[l]] == satdeg[i] ) break;
								}
								if(next[i] == r) {
									aux=satdegheap[nmemb-1];
									satdegheap[nmemb-1]=satdegheap[l];
									satdegheap[l]=aux;
									nmemb--;
									if (satdeg[satdegheap[l]] > satdeg[satdegheap[l/2]]) heapup(satdegheap, l, sizeof(int), max_satdeg);
									else heapdown(l, &satdegheap[nmemb], satdegheap, nmemb, sizeof(int), max_satdeg);
								}
								else satdegheap[l] = next[i];
							} else next[prev[i]] = next[i];
							prev[next[i]] = prev[i];

							// calculate the new sat degree

							satdeg[i]+= (w->weight(i,j) - colorContr[i][c]);
							colorContr[i][c] = w->weight(i,j);
							// insert i in the list with new sat deg

							found = 0;
							for(l=0;l<nmemb;l++) {
								if(satdeg[satdegheap[l]] == satdeg[i]) {
									next[i] = satdegheap[l];
									prev[i] = r;
									prev[satdegheap[l]] = i;
									satdegheap[l] = i;
									found = 1;
									break;
								}
							}

							if(!found) {
								heapoffer(&i, satdegheap, nmemb, sizeof(int), max_satdeg);
								next[i] = r;
								prev[i] = r;
							}

						}
					}
				}
	}

	for (i = 0; i < r; i++)
		delete colorContr[i];
	delete colorContr;

	return nclique;
}

template <>
int WGraphs<int>::dsaturCoverSort(Graph * const g, Weight<int> * w, const int k, int * rank) {
	return dsaturCoverSort_impl<int>(g, w, k, rank);
}

template <>
int WGraphs<double>::dsaturCoverSort(Graph * const g, Weight<double> * w, const int k, int * rank) {
	return dsaturCoverSort_impl<double>(g,w, k, rank);
}




///////////////////////////


class DSaturWeight : public Weight<int> {
	Graph * g;
public:
	DSaturWeight(Graph * g) : g(g) {};
	~DSaturWeight() {};

	int	weight(const int i, const int j);
};

inline int DSaturWeight::weight(const int i, const int j) {
	return g->vertex(i) == g->vertex(j) ? 0 : 1;
}

template <class W>
static inline int dsaturFracColoringSort_impl(Graph * const g, Weight<W> * w, int * r, int ** R, const int k, int ** cover) {
	int returned, i, j;
	unsigned long int ui,mmc_ui;

	mpq_t frac;
	mpz_t den, mmc;

	mpz_init_set_ui (den, 1);
	mpz_init_set_ui (mmc, 1);

	// Calculating the minimun integer 'p' such that p*x[i] is integer for i=1...g->nverts
	for (i = 0; i < g->nverts(); i++) {
		if (w->weight(i, i) > 0) {  // If 'i' is a representative
			mpq_init (frac);
			mpq_set_d (frac, w->weight(i, i));

			printf("x[%d][%d]  converted to %6.4lf", i, i, mpq_get_d(frac));

			mpq_canonicalize (frac);
			mpz_set (den, mpq_denref (frac));
			printf("  den=%6.4lf\n", mpz_get_d (den));
			mpz_lcm (mmc, mmc, den);
		}
	}

	mpz_t cmmc;
	mpz_init (cmmc);

	mpz_set_ui (den, 10);

	mpz_pow_ui (den, den, mpz_sizeinbase (mmc, 10)-2);
	mpz_cdiv_q (cmmc, mmc, den);
	mmc_ui = mpz_get_ui (cmmc); // Common denominator, unsigned long integer */
	printf("mmc=%6.4lf needs %d digits, rounded to %lu\n", mpz_get_d (mmc), mpz_sizeinbase (mmc, 10), mmc_ui);
	fflush(NULL);

	// Allocating memory
	int * numColorAux = new int[g->nverts()];
	memset(numColorAux, 0, g->nverts()*sizeof(int));

	mpq_t mmc_q;
	mpq_init (mmc_q);
	mpq_set_z (mmc_q, mmc);

	mpq_t xi_q;
	mpq_init (xi_q);

	mpz_t xi_z;
	mpz_init (xi_z);

	(*r)=0;
	double xi;
	for (i = 0; i < g->nverts(); i++) {
		xi = 0.0;
		for (j = 0; j < g->nverts(); j++)
			if (w->weight(j, i) > 1.0e-4)  // If 'j' represents i
				xi += w->weight(j, i);

		mpq_set_d (xi_q, xi);
		mpq_mul (xi_q, xi_q, mmc_q);
		printf("x[%d]=%6.4lf   x[%d]*mmc=%6.4lf", i, xi, i,mpq_get_d(xi_q));

		mpz_set_q (xi_z, xi_q);
		mpz_cdiv_q (xi_z, xi_z, den);
		numColorAux[i] += mpz_get_ui (xi_z); // Common denominator, unsigned long integer */
		printf("  numColor=%d\n",numColorAux[i]);
		(*r) += numColorAux[i]; // Counting total colors
	}
	printf("Extended number of vertices: %d\n", *r);

	// Allocating memory
	*R = new int[*r];
	*cover = new int[*r];

	int vertex = 0;
	// Creating R which contains the vertices and its replicas in an alternating order
	for (i = 0; i < (*r); i++) {
		if (numColorAux[vertex] > 0) {
			(*R)[i] = vertex;
			numColorAux[vertex]--;
		}
		else
			i--;
		if (vertex == g->nverts()-1)
			vertex = 0;
		else
			vertex++;
	}
	delete numColorAux;

	Graph * gc = g->complement();
	for (i = 0; i < gc->nverts(); i++) {
		gc->delEdge(i,i);
	}
	Graph * gh = gc->homomorphism(*r, *R);
	DSaturWeight * wh = new DSaturWeight(gh);
	returned=WGraphs<int>::dsaturCoverSort(gh,wh,k,*cover);
	for (i = 0; i < g->nverts(); i++)
		g->delEdge(i,i);

	delete gh;
	delete gc;
	delete wh;

	return returned;
}

template <>
int WGraphs<int>::dsaturFracColoringSort(Graph * const g, Weight<int> * x, int * r, int ** R, const int k, int ** cover) {
	return dsaturFracColoringSort_impl<int>(g,x,r,R,k,cover);
}

template <>
int WGraphs<double>::dsaturFracColoringSort(Graph * const g, Weight<double> * x, int * r, int ** R, const int k, int ** cover) {
	return dsaturFracColoringSort_impl<double>(g,x,r,R,k,cover);
}
